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Abstract. We present a technique for analyzing the full three-dimensional density 
^T^ ' profiles of a planar crystal-fluid interface in terms of density modes. These density 

j/3 ' modes can also be related to crystallinity order parameter profiles which are used 

in coarse-grained, phase field type models of the statics and dynamics of crystal- 
fluid interfaces and are an alternative to crystallinity order parameters extracted from 
simulations using local crystallinity criteria. We illustrate our results for the hard 
'"^ I sphere system using finely-resolved, three-dimensional density profiles from density 

G ' functional theory of fundamental measure type. 
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1. Introduction 

The study of interfaces between the thermodynamically stable bulk phases continues to 
be of great practical and theoretical interest. One main motivation is simply the fact 
that every finite sample of material in some environment contains interfaces. If the bulk 
sample is sufficiently small, the physical properties of these interfaces may even govern 
the behavior of the whole sample. Additionally, surface-specific processes are widely 
investigated in materials science. 

The liquid-vapor interface appears to be the simplest of these interfaces. For a 
simple liquid, there is just one scalar order parameter, the one-particle density, which 
continuously changes when going from the vapor to the liquid. The basic statistical 
mechanical theory of the liquid-vapor interface has already been developed a while ago 
[1] but, of course, there are still open problems such as the question of the correct 
formulation of effective interface Hamiltonians and the behavior of the wavevector- 
dependent surface tension [2], |3] . 

The solid-liquid and solid-vapor interfaces appear to be more difficult to describe. 
Upon solidifying, the continuous translational symmetry of the liquid (vapor) is broken 
to a discrete symmetry describing the crystal. Connected to this, besides density a 
number of additional order parameters are needed which completely characterize the 
bulk crystal and which should vanish when going over to the bulk liquid (vapor) phase. 
Phenomenological attempts to capture the liquid-solid transition from a point of view 
of statistical mechanics have restricted the choice of order parameters to density and 
one crystallinity parameter (models of phase field type |1]) and employed ad-hoc free 
energies of Landau-Ginzburg type to incorporate the transition. A more fundamental 
approach is classical density functional theory (DFT) [1] in which all equilibrium 
properties are governed by a unique free energy functional of the one-particle density. 
In this approach, the crystalline order parameters can be related to the Fourier modes 
of the inhomogeneous density profile. Depending on the level of approximation used, 
only a few modes are needed (as in phase field crystal models |5]) or many (as in 
Ramakrishnan-Yussouff models, weighted density functionals, ...[6] ). In practical 
calculations, however, the density mode expansion is almost always severely truncated. 

From the DFT point of view, the full three-dimensional density profile of a 
solid-liquid interface contains all necessary information about the crystalline order 
parameter profiles. As remarked, with very few exceptions existing DFT calculations 
used restricted parametrizations and none of them attempted a systematic analysis of 
the three-dimensional density profile in terms of modes. In this paper, we present such a 
mode analysis of the full density profile which is general and the resulting modes are also 
quantities which are (in principle) observable in scattering experiments or real-space 
experiments (confocal microscopy). We exemplify the technique using very accurate 
DFT data (using fundamental measure functionals) for the hard sphere crystal-liquid 
interface. It appears also to be interesting to analyze simulation data involving solid- 
liquid interfaces using density modes. Quite often, effective order parameter profiles are 
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measured in simulations through bond-order variables which are, however, difficult to 
capture in theory as they are connected to n-point correlation functions where n is the 
number of next neighbours of a particle or a multiple thereof. 

The paper is structured as follows: In Sec. [2] we introduce the mode analysis in 
general terms. Sec. |3] contains results of the mode analysis of our DFT data. Here 
we briefly introduce the fundamental measure functionals, review the thermodynamic 
observables for the bulk crystal and the crystal-liquid interfaces and then discuss the 
density as well as the free energy modes. Sec. H] contains a discussion on the meaning 
and measurability of the modes and some conclusions. 

2. Mode analysis 

Consider the density distribution pcr(i" = {x, y, z)) in a perfect crystal. It can be Fourier 
expanded into a discrete sum, 

Pcr(r) = ^Pj exp(iKj • r) , (1) 

j 
where Kj denotes the set of all reciprocal lattice vectors (RLV) and Pj are the associated 
Fourier amplitudes. Owing to the symmetry of the crystal under consideration, not all 
Pj are independent of each other. It is convenient to group the Kj in shells with index m 
where all Kj belonging to one shell can be transformed into each other under the discrete 
symmetry group of the crystal under consideration (and thus all Pj associated with these 
Kj are equal). As an example, for the face-centered cubic (fee) crystal the reciprocal 
lattice is of body-centered cubic {bee) symmetry. Let a be the side-length of the cubic 
unit cell of fee and correspondingly b = 271 /a the side length of the cubic unit cell of bcc 
in reciprocal space. The reciprocal basis is given (in Cartesian coordinates where the 
axes span the cubic unit cell in reciprocal space) by Bi = 6(1, 1, —1), B2 = 6(1, —1, 1) 
and B3 = 6(— 1, 1, 1). An arbitrary RLV is a linear combination of the B,. The shells 
are characterized by a triple (m, n, k) of natural numbers and the Kj belonging to this 
shell have Cartesian components 6(±?7i, ±n, ±fc) and permutations thereof. Thus, if 
m,n,k are mutually distinct, there is a maximum of 48 RLV in one shell. The shells 
with lowest modulus are given by (1,1,1), (2,0,0) and (2,2,0). A hsting of the RLV 
triples up to shell 15 is given in Ref. [7] (Table I). 

We prepare a planar crystal-fluid interface with interface normal in z-direction 
and density profile p{x,y,z). It is tempting to parameterize the density profile using 
the RLV as before, 

p(x, y,z) = Y^ exp{iKj ■ r) pj{z) , (2) 

j 
but now with a z-dependent Fourier amplitude Pj{z) for the reciprocal modes. We call 
them (density) modes or order parameter profiles. One expects that upon crossing the 
interface coming from the crystal side, all Pj{z) go to zero (for nonzero Kj). Only for 
Kj = the value for the associated mode crosses from the average density of the crystal 
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Figure 1. (color online) Hard sphere crystal-fluid interface in [100] orientation, (a) 
Laterally averaged density profile, (b) Fourier transform of laterally averaged density 
profile. Clearly visible are the "lumps" around multiples of the wave numbers 47r/a 
(corresponding to the periodicity a/2 in the density profile). The back transforms of 
the individual lumps define the mode profiles. 



to the average density of the fluid. As discussed before, in the bulk crystal all RLV 
within one shell with index m are degenerate. At the interface the degeneracy of the 
RLV in one shell is lifted, and we introduce an index n which distinguishes the possible 
values of the z-component of the RLV. Thus the decomposition becomes 






P(a;, 2/, z) = 2_^2^Pmn{z) exp (i(Kj) 



(3) 



The sum over j is now only for those RLV within shell m which have common value of z- 
component, as expressed by the index n. In the literature, such a decomposition has been 
used to parameterize the full 3d density profile p{x, y, z) using only the leading mode 
Pinilin] in order to facilitate a simplified order parameter description of the crystal-fluid 
interface. Consistency with the bulk crystal solution pcr requires Pmn{.z — )■ — cxd) = P^, 
thus in a leading-mode description either the crystal description is also truncated [10] 
or higher-than-leading modes are parametrized proportional to the leading mode [9]. 

For a given 3d density profile p{x,y,z) the extraction of the modes pinni.^) does 
not appear to be straightforward since they can not be projected out using a Fourier 
integral. Consider the Fourier transform of the 3d density profile: 



p\kxi ky, kz 



LxLy 2_^ 2^ 



5^ (k± + (K_Lj)mn)Pmn {kz + {Kz 



(4) 



Here, the lateral dimensions of the system with interface are denoted by L^ and Ly. 
The wavevector in the interface plane, /c_l = {kx, ky), is restricted to the discrete values 
allowed by the RLV, (K±j)mn- The Fourier transform of the modes, Pmn, are shifted 
by the z-component of the associated RLV, {Kz)mn, which is by our convention fully 
specified by the indices m,n. Hence we see that the Pmn{kz) can be viewed as being 
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centered around —{Kz)mn and are well separable from each other as long as their range 
in Fourier space does not extend further than the separation between two neighbouring 
{Kz)mn- This in turn is given by 27r/c?[o] with d[o] being the z-separation of lattice planes 
in the orientation [o] . This condition on the mode separability is a very physical one: 
it means that the density modes should vary across the interface much more smoothly 
than the actual density profile p{x, y, z). 

We illustrate the mode separability using explicit DFT results for the crystal-fluid 
interface in the hard sphere system, for more details see Sec. |3] below. For the [100] 
orientation, Fig. [1] (a) shows the laterally averaged density profile 

\ pLx rLy 

P-^Az) = -f—j- j dx dyp{x,y,z), (5) 

J^x-l^y Jo Jo 

and Fig.[T](b) its Fourier transform, corresponding to p(0, 0, kz). This Fourier transform 
consists of localized "lumps" around the values i-Air/a {i = 0, 1, 2, . . .) which mark lattice 
planes separated by the distance d = a/2 (which is what we expect for the laterally 
averaged density of an fee crystal). The inverse Fourier transform of the isolated "lumps" 
(centered around 0) gives the associated modes in real space. As indicated graphically 
in Fig. [1] (b), the Fourier transform of a particular mode can be isolated by cutting out 
the Fourier transform p(0, 0, k^) symmetrically around i ■ 47r/a with a width An/a. We 
formalize this a bit more generally as 

Pmn{kz) = p{Kx, Ky, kz + {Kz)mn) g{kz) (6) 

where (K^ j)mn = {K^, Ky) denotes the in-plane components of the RLV of interest for 
mode {mn) (again, there is an additional degeneracy, expressed by j). The function 
g{kz) is a suitable cutoff function, varying monotonically from 1 at around fc^ = to 
zero near kz = 27c /d[o] and with the additional properties 

27r' 



9{kz) = \kz\ > , 

9{kz)-g{-kz) =0, (7) 

2n' 



9ikz)+g[ -^-kz] =1 [kzE 

d\o] 



"fol 



Note that these conditions are required to have (EI) consistent with (j4j). A practical 
choice (which violates the first condition a tiny bit) is 



g{kz > 0) = - ( 1 - erf ( a 



kzd[o] _ ' 



(8) 



7r 

where a tunes the width of the error function kink. Note that for a — > cxd the cutoff 
function becomes g{kz > 0) = 6{n/dio] — kz), literally corresponding to the graphical 
procedure indicated in Fig. [1] (b). The disadvantage hereby is that through a hard cutoff 
in Fourier space unwanted oscillations arise in real space which simply might be due 
to insufficient resolution of the density profile in real space and are not due to physical 
oscillations of the modes. We found a > 3 a convenient choice. 



Mode expansion for the density profile of crystal-fluid interfaces: Hard spheres as a test case 6 

We had introduced the modes Pmn{z) as belonging to the RLV shell m but differing 
in the z-component of the associated RLV, as indicated by n. If that 2;-component is 
zero, the mode will be purely real. If for a particular n that z-component is nonzero, 
the mode will be in general complex since from Eq. ([6]) it can not be expected that the 
Fourier transform of the mode is symmetric in kz- The reality of p{x, y, z) requires that 
in the same shell m there is a mode Pmn{z) for which the 2;-component of the associated 
RLV is minus the previous one and which is the complex conjugate oi Pmn{z). Hence 
can define real-valued modes through separating real and imaginary part, 

-(_ / N Pmn\^ ) ' Pmn\^ ) __ / \ Pmn\^ ) Pmn\^ ) 

Pmn\^) 7, ' Pmn\^) 



2i 

and the Pmn(^) have the obvious interpretation of phase shifts of the associated density 
mode oscillations across the interface. 

3. Results for the hard sphere crystal-fluid interface 

3.1. Fundamental measure DFT 

We apply the mode expansion technique introduced in the previous section to DFT 
results for the full crystal-fluid interface density profile. DFT is built on the existence 
of a unique free energy functional of the one-particle density field p(r), 

j^[p] = r^ip] + j^^^ip] , (9) 

(3r'^[p]= /"rfVp(r)(ln(p(r)A3)-l) (10) 

which can be split into the exactly known ideal gas part J-""^ (A is the de-Broglie 
wavelength, /3 = l/ksT is the inverse temperature) and a generally unknown excess 
part J-"*^^. The equilibrium density Peq(i") in the presence of an external (one-particle) 
potential V^^^lr) is then given by 

6J^[p] 



6piv) 



= p-V^-\r), (11) 

P=Poq 

where p is the imposed chemical potential (e.g. by requiring a certain bulk density 
far away from the region where the external potential acts). For the description of 
the equilibrium density profile between crystal and fluid, V^^^ = 0. Thus the bulk 
crystal far away from the interface on one side appears as a self-sustained, periodically 
inhomogeneous fluid at the coexistence chemical potential /icoex- 

For the hard sphere system, fundamental measure theory (FMT) allows the 
construction of very precise functionals [HI |12l Il3]. Essentially, FMT postulates an 
excess free energy with a local free energy density in a set of weighted densities n^: 

^-[p]=/3-i [ dh^n^iv)). (12) 



The weighted densities are constructed as convolutions of the density with weight 
functions, nQ,(r) = p * w°'{r). The weight functions reflect the geometric properties 
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of the hard spheres. For one species, the weight functions include four scalar functions 
w° . . . w^, two vector functions w", w^ and a tensor function w* defined as 

in 111 

w^ = e{R-\v\), w'' = 5{R-\r\) , w^ = ^ , w" - 



AttR ' 47ri?2 



,2 r 1 w 



2 



w' = — 5(i?-|r|) , w 



rl ' ' '' ' 4:7iR 



wl, = ^6{R-\r\). (13) 

Here, R = cr/2 is the hard sphere radius. Using these weight functions, corresponding 
scalar weighted densities uq-.-u^, vector weighted densities ni,n2 and one tensor 
weighted density Ut are defined. In constructing the free energy density $, arguments 
concerning the correlations in the bulk fluid, certain geometric consistencies and 
arguments pertaining to strongly inhomogeneous systems are used [T3] : 

H{n[pir)]}) = -no ln(l - ^3) + ^iK) !^l!!lZJli:J^ + 

1-723 



3 (-722 n2 • n2 + n2,int,ijn2j + ^2 rit^ijUtji - nt^ijntjknt,ki) /. ,. 

<^2 ri3 ;' ' — ^^^-^ -^^ — — ■ 14 

lo7r(l — n^j^ 

Here, ipi{ns) and v'2(^3) are functions of the local packing density ri3(r). The most 
accurate functional to date is the tensor version of the recently introduced White Bear 
II (WBII) functional [H] and corresponds to the choice 

,^.^-l \ ^^3 - nj + 2(1 - ns) ln(l - ^3) 

3^3 

2^3 - 3^2 + 2nl + 2(1 - ^3)2 ln(l - ^3) 

"' = ' H ■ 

This functional has proved to be very reliable for strongly inhomogeneous fluids [151 [IS] 
and also wall-fluid surface tension data at densities close to freezing excellently compare 
with recent simulation data [T7] . 

3.2. Thermodynamic properties of the crystal and the crystal-fluid interface 

In equilibrium, hard spheres in their crystalline state form an fee lattice, although the 
free energy difference to hep is very small, about 10~^ kBT per particle [18]. Bulk 
coexistence properties (densities, chemical potential, pressure, free energy) have been 
widely investigated by simulations and there seems to be consensus on the respective 
values with an error margin below 1%. The bulk crystal has been investigated in 
Ref. [21] using full minimization of the FMT functional discussed above and very good 
agreement has been found with simulations. The coexistence density on the liquid 
side is overestimated by 0.5%, this results in an overestimation of the coexistence 
chemical potential and pressure of about 2%. The relative vacancy concentration 
nvac = 2 ■ 10~^ at coexistence is too small by a factor of 10 compared to simulations, 
however, an assessment of possible systematic errors in the simulations is not available. 
The FMT results for n^ac improve significantly on results from Ramakrishnan-Yussouff 
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1 From Ref. [19]. 

2 Free energy for PcrO'^ = 1.041 [I2j using an improved fit in the form of the Speedy 
equation of state from Ref. pT] . 

3 Chemical potential for p^a^ = 0.940 [19j from the Carnahan-Starling equation of 
state. 

^ From Ref. EOl. 
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From Ref. 



Table 1. Thermodynamic properties of hard sphere crystals and crystal-fluid interface 
from FMT in comparison to simulation data. 



and weighted density functional theories where rZvac ~ 0.1. Also, in contrast to these 
earlier approaches, the FMT crystal density profiles agree with simulations with respect 
to the width of the crystal peak and the anisotropy in different lattice directions |2T] . 
See Table [1] for a summary. 

The determination of crystal-fluid surface tensions in simulations has not yet 
resulted in consistent data across the literature. As an example, for the surface tension 
in [110] direction integration methods give /37[ioo]0"^ = 0.58 (cleaving, Ref. [22]) and 
0.64 (umbrella sampling, Ref. p3])- Through the capillary wave method, interfacial 
stiffnesses are accessible but the conversion to tensions is not without difficulties owing 
to the finite cutoff in the stiffness and tension expansion in cubic harmonics j24j . Here, 
Laird et al. find /97[ioo]C"^ = 0.58 (expansion up to order two, Ref. [25]) and newest 
results give 0.64 (expansion up to order three, Ref. [21]). The FMT results (also from 
Ref. [21]) are larger by about 7% and the anisotropies compare well to the simulation 
results, see Table [1] It is perhaps not so surprising that the surface tension in FMT is 
higher than in simulations since the FMT functional does not capture very long-ranged 
fiuctuations present in capillary waves. Note that previous theories based on (i) severe 
approximations to the free energy functional and {ii) very restricted minimizations in 
density space have produced similar results [261 122], consequently the agreement in 
numbers has to be regarded as fortuitious. See also below for the discussion of the 
density modes which puts simple density parameterizations in perspective. 
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Figure 2. (color online) Cuboid unit cells for different orientations of the interface 
normal (z-direction). Different layers with z — const, are given in different colours. 
For [210], layers n = 7... 11 are omitted for clarity. They can be obtained from 
layer k = 12 — n using the reflection of coordinates x — )■ v^a — x for all particles 
in the layer. The unit cell for the [211] orientation is obtained from [111] using the 
axis relabeling a; — > j/, y — ^ z and z — >■ a;. The z-separation between lattice planes, 



needed in Eqs. 1^ and (O, is given by d[ioo] = Q.z/2 = a/2, d[ioo] 
d[iii] = a^/3 = a/V3, dpn] = a^/G = a/{2\/6) and d^io] = a^/lO 



: a,/2 = a/{2V2), 
a/(2V5). 
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Table 2. The assignment of reciprocal lattice vectors to the main shell index 
m — ... 3 and z-component index n for the [100] interface. Note that modes with 
index n and n are related by complex conjugation. 



3.3. Density modes 

We analyze the density profiles of crystal-fluid interfaces in five different orientations of 
the interface normal: [100], [HO], [HI], [211] and [210]. The density is given for a cuboid 
with edge lengths L^,, Ly and L^ which contains the fluid in the middle {z ~ L^jl) and 
the crystal phase at the boundaries (z ~ and z ~ L^). Lx[y\ are given by the edge 
lengths in 2:[y]-direction, 0^;^^], of the smallest cuboid unit cell of the crystal which has 
the desired orientation in 2;-direction. Typically we chose L^ = 320^. The crystal cuboid 
unit cells are depicted in Fig. O Whereas they are still comparatively easy to imagine 
for the [100], [110], [HI] and [211] orientations, it is less trivial for the [210] orientation. 
As explained in Sec. O the density modes Pmn{z) are labeled by their RLV shell 
index m and by n for the possibly different 2;-components of the RLV. The specific 
RLV assignment within one shell for the index n then depends on the chosen interface 
orientation. As an example, we show the assignment for the [100] interface in Table El 
We call the mode ^11(2;) the leading crystallinity mode since the associated RLV is in 
the lowest shell and in the bulk crystal the associated Fourier amplitude Pi is largest. 
The meaning of pu (z) is linked to the strength of density oscillations of close-packed 
planes across the interface into the bulk liquid. For the [100] interface, the normal to 
the close-packed planes is not in z-direction. Another mode of significance is the mode 
P22{z) which is linked to the strength of density oscillations of square arrays of particles 
which are in the x-y-plane. It is the leading Fourier component of the lateral density 
average (often shown in simulation works) when the average density mode poo{z) has 
been subtracted. 
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Figure 3. (color online) Density mode profiles Pmn{z) in units of a^^ up to m = 4 
for the [100] interface, (a) Real part of profiles, note that the density mode poq{z) has 
been rescaled and shifted, see tick labels at right y-axis. (b) Imaginary part of profiles. 
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In Fig. [3] we show the modes up to m = 4 for the [100] interface. In Fig. |3] (a), 
the real part of the mode profiles is shown, including the rescaled density mode ^00(2;) 
(fat dots). The density mode shows a depletion zone in front of the crystalline phase. 
In previous work on the crystal growth at walls [28], a depletion zone has been found 
in instantaneous profiles and attributed to the finite velocity with which the interface 
moves. The absence of the depletion zone for equilibrium profiles as seen in Ref. [28] 
can be attributed to the averaging procedure: there the density was averaged between 
minima of the laterally averaged density profile. Due to the phase shift of the density 
oscillations across the interface, such an averaging procedure does not give our density 
mode poo{z); we verified that using the procedure of Ref. [28], also the depletion zone 
vanishes. Further inspection of Fig. [3] (a) shows that the leading crystallinity mode 
Pu{z) shows a smooth, monotonic behavior which corresponds to the expectations on 
a phase-change order parameter profile and which could be used to define an interface 
location. However, this obviously involves some arbitrariness: the kinks in higher modes 
m > 2 are clearly positioned deeper in the crystal phase and the mode P22{z) (leading 
crystallinity mode of the lateral density average) differs also in the qualitative shape 
by showing a local minimum (similar to the depletion zone in poo{z)). Furthermode 
we note that the difference between the interface location as extracted from the leading 
crystallinity mode ^11(2;) and the density mode poo{z) differ by about one unit cell length 
a. We discuss this further below. - In Fig. [3] (b) we show the imaginary part of the mode 
profiles, corresponding to the phase shift of oscillations across the interface. Positive 
values indicate that the distance between maxima becomes larger in the interface region, 
see also Fig. [1] (a) where this is clearly visible for the laterally averaged density profile. 

As it is the case for the surface tensions, also the modes are not rotationally 
invariant. For a new orientation, one takes the rotation matrix which transforms the 
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Table 3. The assignment of reciprocal lattice vectors to the mode indices m,n 
for the first two main shells and for the orientations used (except [100]). ifs[y,2] are 
given in units of 2-K/anc\^y^z] where flajf^.z] are the side lengths of the cuboid unit cells of 
Fig. [2] with the desired interface orientation in z~direction. Within a shell m — const., 
the RLV are sorted with increasing |i^z|, for simplicity n and n have been identified. 
Marked with an asterisk (*) are the monotonously varying leading modes shown in 
Fig. [5] below. 



coordinate system of the [100] unit cell to the cuboid unit cell with the desired orientation 
(see Fig. |2]) and applies it to the RLV. The RLV within a shell m need to be regrouped 
with respect to their new 2;-components (see Table [3]). Whenever a z-degeneracy is 
lifted when going to another orientation (such as for the leading mode m = 1 which is 
completely degenerate in [100] orientation), there will be new modes. As an example, 
in Fig. Hlwe show the modes up to m = 2, both real and imaginary part, for the [111] 
interface. The shell m = lis split into two modes where the associated RLV in reduced 
units are given by (0, 0, 3) (for pu^z)) and (0, 2, 1) (for pi^^z)). (The reduced units are 
defined by normalization of RLV component x[j/, z] to the value 2T[/ax[y,z] where a^iy^z] 
are the side lengths of the cuboid unit cell, see Fig. El) As can be seen in Fig. HI mode 
Pii{z) is similar in shape to the same mode for [100] and corresponds to the decay of 
the density oscillations of the close-packed x-y-planes. However, the decay of density 
oscillations of the other close-packed planes which are not perpendicular to the interface 
normal is different, as the behavior of pi2(^) shows. 

We also found that for the other orientations there is always one leading mode 
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Figure 4. (color online) Density mode profiles Pmn{z) in units of a^^ up to to = 2 for 
the [111] interface. Note that the density mode pm{z) has been rescaled and shifted, 
see ticks at right y-axis. 



{m = 1) which shows the monotonous decay of the leading mode for the [100] orientation, 
for the associated RLV see Table [31 Also the shape is roughly similar for all orientations, 
and can be roughly fitted to a tanh-profile with width a (see Fig. [5] (a)). Also the 
average density mode qualitatively has the same shape (especially the depletion zone) 
for all orientations, see Fig. [5] (b). Therefore our previous finding holds regarding the 
clear separation between interface position, as determined from the leading crystallinity 
mode, and the interface position, as determined from the average density mode. This 
separation is approximately of size a ~ 1.6 a which is not small. So, coming from 
the liquid side, the hard sphere fluid first orders and then densities, in contrast to the 
time sequence in homogeneous nucleation where nuclei form by first densifying and then 
ordering [29] . 

Since the interparticle interaction is short-ranged, one would expect an (oscillatory) 
exponential decay of the modes into the respective bulk phases. An analysis within 
Landau theory predicts an upper bound of the decay length of mode Pmn{z) of 
l/|(K_Lj)mn| [3Q]. A more refined analysis taking into account the direct correlation 
function of the bulk liquid [H [31] also finds that the decay lengths decrease with 
increasing modulus 1/| (K_Lj)mn| but the magnitude is enhanced by a factor of 2. . . 3. We 
defer a detailed analysis of the decay lengths to future work and merely remark that the 
aforementioned simple rules on the decay lengths do not seem to capture the qualitative 
importance of modes in the approach to the bulk phases, see Fig. [3l (For example, poo 
and P22 should have the same decay length but the corresponding amplitudes are very 
different in the approach to the liquid phase. Furthermore, the different ^-position of 
the mode kinks is very important for the relative importance of the modes at a given 
location on the z-axis.) 

In view of all the features we have seen in Figs. [3]-[5|and discussed above, it is clear 
that restricted density parametrizations can not be very reliable. This has been shown 
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Figure 5. (color online) (a) Leading mode (in units of cr^'^ for all five orientations and 
a fit of tlie form pfit{z) — 0.45(1 + tanh[(z — zo)/(wa)]) with z^/a — 6.35 and width 
w = 1.0. (b) Average density mode for all 5 orientations, for comparison the fitted 
leading mode, rescaled and shifted to match the asymptotic densities, is given. 
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already for the hard sphere crystal-hquid interface evaluated with a simpler functional 
of weighted density type. Full minimization [32] reduces the surface tension by a factor 
of 2 compared to a parameterization [26]. Also for FMT functionals quite elaborate 
parameterizations have been tested in Ref . [33] . The average surface tension is still larger 
by about 30% than the results given in Table [H and the ordering 7[ioo] > 7[iio] > 7[iii] 
is not reproduced. 



3.4- Free energy modes 

The primary focus of the analysis in the last section was on the density modes which 
are observable quantities (see also Sec. H] below). Nevertheless, the mode analysis can 
also be performed for the free energy density / or the grand free energy density u using 
the DFT results. The particular significance of the leading mode ujqq{z) of the grand 
free energy density is the connection to the surface tension: 



7 



dz{ujQQ{z) + Pc 



(16) 



In Fig. |6]we show the leading mode both for the free energy density (a) and for the grand 
free energy density (b). The leading free energy modes for the five different orientations 
actually look very similar and almost monotonously connect from the higher free energy 
density of the solid to the lower free energy density of the liquid (save for a small 
hump on the crystalline side). It becomes clearer in the plot for the grand free energy 
density that there is indeed a broadening of the interface (in free energy terms) in the 
sequence [100]-[211]-[110]-[111]-[210]. This is also in rough agreement with the width of 
the leading crystallinity modes shown in Fig. |5] (a). 
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Figure 6. (color online) Leading free energy modes, (a) Free energy density, (b) 
Grand free energy density, shifted by the coexistence pressure. 



Modes other than the leading one appear to be without much physical content. 
However, these may be used to gauge more phenomenological approaches such as phase 
field crystal models [5]. 



4. Summary and Conclusion 

In this paper, we have presented an analysis of the full 3d density profile of a crystal- 
liquid interface in terms of density modes associated with the reciprocal lattice vectors 
of the bulk crystal. We have exemplified the extraction of the density modes on finely- 
meshed density functional data obtained for the hard sphere system and found a number 
of interesting results: 

• a separation of about 1.6 cr between the interface location as determined by the 
average density and the interface location as determined by the leading crystallinity 
mode 

• a density depletion zone just in front of the bulk crystal 

• strongly non-monotonous mode profiles also for next-to-leading modes 

We suggest that simulation data on crystal-fiuid interfaces should be analyzed for these 
modes. According to Eq. IQ, the Fourier transform Pmnikz) of the mode {mn) is 
obtained by "cutting out" the fc^-component of the Fourier transform of the 3d density 



profile p{kx, ky, kz) around the reciprocal lattice vector k 



K 



] )mn 



(different j for 



an index pair (mn) signal mode degeneracy). Practically, this means just averaging 
(exp(i(Kj_j)mn ■ i"±)) in the x-y-plane as a function of z. The average would guarantee 
that the statistics will be good enough for standard runs that one obtains smooth 
Fourier transform and can perform the extraction procedure. Alternatively, also the 
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Fourier transform in z-direction can be sampled directly. This procedure appears to be 
less prone to statistical uncertainty than sampling a full 3d profile. 

In principle, the Fourier transforms of the modes pmn{kz) are observable quantities 
in scattering experiments. Practically, one would need to irradiate a volume containing 
crystal, liquid and the interface, fix the lateral momentum transfer to a particular 
(K^ ,,)m„ and finely scan the momentum transfer in z-direction. In view of additional 
effects like incoherent background and finiteness of the volume with a well-defined 
interface orientation, this seems to be a very demanding task. Previous studies on 
surface melting using LEED [M] for instance only obtained some integral information 
on the modes by recording an smaller Debye- Waller factor (enhanced decrease in peak 
intensity) for a certain reciprocal lattice vector due to the decay of the associated mode. 
The theoretical analysis done at that time was very rough and could be improved using 
the techniques presented here. 

It is perhaps easier to use confocal microscopy for the observation of colloidal 
crystals [2H1 EHl |36] and to do a planar averaging as described for the analysis of 
simulation data. In fact, confocal studies on crystals have previously frequently 
employed techniques which are borrowed from simulation such as the bond-order 
analysis. As remarked in the Introduction, observables connected to bond order 
are difficult to handle theoretically, so the mode analysis of simulation and confocal 
experiments will be beneficial for theory development. 
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